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Abstract 

A numerical method close to the Strutinsky procedure (but better) is proposed to calculate the deformation energy of nuclei. 
Quadrupole (triaxial) deformations are considered. Theoretical as well as practical aspects of the method are reviewed in this 
paper. A complete fortran program illustrates the feasibility of the method. Thus, this code will constitute a useful "ready 
tool" for those which deal with numerical methods in theoretical nuclear physics. 
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I. INDRODUCTION 



There are two methods which allow to determine the equilibrium shape ( ground state) of the nuclei: The constrained 
Hartree-Fock method and the so-called macroscopic-microscopic method. Though the latest generation of computers 
is able to perform very complicated calculations, in terms of running time, it is no so obvious to make systematic 
calculations for a large number of nuclei. A good alternative is to use the Strutinsky method. The latter consists of 
associating the classical liquid drop model with some shell and pairing corrections built from a realistic microscopic 
model. Based on such a model, we present a numerical method with its associated fortran program . The potential 
energy of deformation is deduced as a function of the shape of the nucleus. Triaxial (quadrupole) shapes are considered 
in this work. The three semi axes of the ellipsoid are in fact connected to the both Bohr parameters which are actually 
used in the calculations. 

The different steps of calculations are: 

i) The energy of deformation of the liquid drop model is first calculated [4]. 

ii) The Schrodinger equation of a microscopic Hamiltonian is built and solved to obtain eigenvalues and eigenvectors. 
In fact we use the Fortran program named "triaxial" already published in cpc. The microscopic model is explained 
in details in this paper and also in Ref. [3] . 

iii) The semiclassical energy is deduced from the same Hamiltonian as (ii) is calculated on the basis of the Wigner- 
Kirkwood expansion [5]. 

iv) The shell correction is deduced as the difference between the sum of single-particle (point (ii)) energies and the 
same quantity smoothed semiclassically (point (iii)). 



II. POTENTIAL ENERGY OF DEFORMATION IN THE LIQUID DROP MODEL 

We use the so-called macroscopic-microscopic method [l[ to evaluate the potential energy of deformation of the 
nucleus. This method is based on the liquid drop model plus shell and pairing corrections deduced from a microscopic 
model 0. 

The deformation (or potential) energy of the nucleus is defined as the difference between the binding energy of the 
deformed drop and the non-deformed drop (nucleus). 

E(P) = E LD (/3) - E LD (0) (1) 

Here f3 is a set of parameters defining the deformation. The case /? = represents the spherical shape (i.e., the 
non-deformed nucleus). We recall that in the liquid drop model, the minimum is always obtained for the spherical 
deformation. This involves E(/3) > 0. Of course, the liquid drop or weizsaker formula model contains several terms, 
but only two depend on the deformation of the nucleus, namely the surface and the coulomb energies. Consequently, 
the other terms do not survive in the difference given by Eq. (|T|) . The liquid drop energy reads [4j 
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with ro = 1.275 fm and e 2 = 1.4399764 MeV . The quantities B s , et B c are the surface and the coulomb contributions. 
It is to be noted that B s and B c are dimentionless and normalized to the unity so that the deformation energy of the 
non deformed nucleus is equal to zero (i.e., E(0) = 0). The reduced fissility has been determined empirically 
C = 52.8(1 - 2.84/ 2 ), 1= (N- Z)/(N + Z) 

For triaxial ellipsoidal shape with semi axes a, 6, c, the coulomb and surface contributions are deduced analyticall 
with the help of elliptic integrals of the first and second kind F(ip, k) et E(tp, k) so that if a ^ b > c, we will have 

B c = F[ip, k)a 2 b 2 c 2 R^/{a 2 - c 2 ) 1/2 , (3) 
sin^ = (1 - cVa 2 ) 1 / 2 , k 2 = (a 2 - b 2 )/{a 2 - c 2 ) 



B s = „ ■ 
2 



\ Rq 2 {c 2 + b(a 2 - c 2 )-V2 [(fl 2 _ C 2 )E{ ^ W) + C 2 F[ ^ k ,)] } (4) 

sin ip = (1 - cVa 2 ) 1 / 2 , k' 2 = a 2 b- 2 {b 2 - c 2 )/(a 2 - c 2 ) 

the condition of the volume conservation of the nucleus being abc — Rq (equal to r^A). with this condition it is clear 
that only two deformation parameters are necessary to specify the shape of the nucleus. The Bohr parameters (/3, 7) 
are more commonly employed in this type of calculation. For moderate deformations, the link between the semi axes 
and the Bohr parameters is given in Ref. The elliptic integrals are evaluated with Gauss quadrature formulae 
with 64 points. 
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III. SHELL CORRECTION 



According to the Strutinsky prescription, the shell correction to the liquid drop model is defined as the difference 
between the sum of the single-particle energies of the occupied states and the " smoothed part" of the same quantity: 



SE = £ e *-E e * 

occupied 

In fact the Strutinsky procedure is done in such a way that the smoothed sum does anymore contains shell effects so 
that the above difference represents only the contribution due to the shell structure. In the Strutinsky's method the 
smoothed sum is derived through the smoothed density of states [l[ 



(£ e = / e 3i\/, 7 (e)* 

\ / strutinsky J 

Here, g M 7 is the level density and M and 7 are respectively the so called order and smearing parameter of the 
Strutinsky's procedure. The major defect of this method is that generally the results are usually more or less 
dependent on these two parameters. A method to diminish this dependence is to use the plateau condition, however 
in the case of finite wells this is not systematically guaranteed. In this respect, it has been demonstrated in ref. 
that the level density given by the Strutinsky method is nothing but an approximation of the semiclassical level 
density, i.e. a quantum level density from which the shell effects have been washed out. Consequently, even though the 
Strutinsky is simpler in practice, it is more interesting to work straightforwardly with the semiclassical density because 
the problem of the dependence on the two above parameters is in this way avoided. Thus, it is simply recommended 
to perform the smoothing procedure with the semiclassical level density. The previous formula becomes in this case: 



eg sc (e)de 



where g sc is the semiclassical level density. 

Even though it is not so obvious to derive a semiclassical density of states from a given quantum Hamiltonian, there 
is for our case a rigorous solution (in the sense where it the same quantum Hamiltonian which is " treated" semiclassi- 
cally). Indeed, for exactly the same Hamiltonian employed to determine the eigenstates, the semiclassical level density 
is deduced following the Wigner-Kirkwood method. The latter is based on the Thomas-Fermi approximation plus a 
few corrections appearing as a power series of In this theory, the particle-number is expressed as a function of 

the Fermi level as follows [5]: 
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where V and S are the central field (including the coulomb potential for the protons) and the spin-orbit field (see 
Ref. H). The classical turning points are defined by A — V(r^ c ) = 0. The domain of integration is defined by: 



D : V(~?) < A = V(r 
The semiclassical level density is thus derived as follows: 

dN(e) 



9sc{e) 



de 
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(7) 



and the semiclassical energy is therefore: 



E £ i) =E SC = eg sc (e)dt (8) 
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As already mentioned, the Fermi level is obtained from the following equation: 

N(X) = N (9) 
where iVo is the particle-number (neutrons or protons). 

The semiclassical energy which is of course free from shell effects can be cast under a power series of (1/fi): 

E sc = XN Q - (E% + E\ + E°) - (E s ° + E?°) (10) 
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where for example E°_ 3 contains the term (1/ft) 3 , etc... Here SO means the term related to the spin-orbit interaction. 
The expressions of E® et Ef° are very complicated and become simple only for the non-deformed case (spherical 
shape) . The importance of these terms decreases rapidly. The reference [j| gives the following percentages with respect 
to the total semiclassical energy: XN ~ 30.5%, E°_ 3 ~ 66%, £° x ~ 2%, E s f ~ 1.25%, E% ~ 0.05%, Ef° ~ 0.2%. 
In addition, it is to be noted that the "active part" due to the deformation is even smaller. For this reason the 
contributions E® et Ef ° (which are not given explicitly here) are simply approached by their values for the spherical 
shape: 
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The different integrals (ITT|) . (fT2")l . (fT3f are calculated by the three dimensional Gauss-Legendre quadrature formulae. 

The set of lattice points must verify Eq. ([5]). In fact, for convenience, in each direction, each interval is divided in 

elementary intervals in which the quadrature formula is applied with a restricted number of nodes. The number of 

points is increased in such a way to obtain stable numerical results. 

The Fermi level is not determined straightforwardly from Eq. (0), but solved as follows: 

From: 



E x , 



eg sc (e)de = / e^J^-de 



a simple integration by parts gives: 



E SC (X) = XN(X) - / N(e)de with g sc (-oo) = 0. 

J — oo 

with the condition of the Fermi level N(X) = No, we will have 
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E SC (X) = XN - / N{e)de 

J — oo 

the differentiation with respect to A gives 
= N - N(X) = 

This means that for the constraint N(X) — Nq, the value of A is the one which makes E SC (X) minimum. Consequently, 
for a fixed Nq it is sufficient to look for this minimum with the help of Eq. (flOl) (this is what is done in the fortran 
program) without employing subsidiary Eq.Q. Knowing A, the correctives terms E® et Ef° are deduced in the 
spherical approximation (as mentioned before, the dependence on the deformation being very small for these terms) . 
Unlike the previous case, the integral (fT4"|) and (fT5|) are one-dimensional and are also treated by Gauss-Legendre 
formula. It is to be noted that the nodes of the quadrature do not make any problem for the term (A — V) , i.e., 
we have always A ^ V(r no d e )- 



IV. DETAILED EXPRESSIONS OF (VS) AND VT 

Expressions (VS) 2 and V 2 V are derived analytically, for (VS) 2 the result is: 

^ 2 = TT^^{i F{a> + i m + i F<c} } (16) 

„ ( y/rj- 1 ^\ x 2 y 2 z 2 x 2 y 2 z 2 

E=exp {-^r^)' r2 = ^ + -^ + ^ ri = ^ + ¥ + ^ (17) 

F{e) = { J= - + X (l 8) 
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with S(r) = k/ [1 + exp(R so L so /a so )] . It is worth to note that the spin-orbit coupling constant k and present 
work) is related to kj of Ref. Q by the following equation: 

= (-2Mc 2 /h 2 c 2 ) k ~ -0.0482k (19) 

This is due to the fact that in these references, the spin-orbit constant is not defined in the same way. For V 2 V we 
have: 
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1 + exp(ity Lv/av) 

E = cxp(R v L v /a v ) (24) 

F(e) being defined by Eq. (TTSJ). Finally, the shell correction is calculated by replacing the Strutinsky's level density 
by the semiclassical energy: This leads to: 

5E= J2 ti-F sc (25) 

i occupes 

The shell corrections are calculated separately for the neutrons and the protons and then added to obtain the total 
shell correction. 
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V. PAIRING CORRECTION 



We have took into account the pairing correction via the simple BCS approximation. The Fermi level A and the 
jap parameter A are solved from the well known system of coupled equations: 

2 ^ 1 



G 




(27) 



£fe - A 



fe=l V v /(e fc -A) 2 + A 2 



In these equations G is the pairing strength and tk the eigenvalues of the microscopic Hamiltonian. The upper index 
N p of the sums represents the number of pairs of quasiparticles actually taken in the calculations (with Np/2 above 
and below the Fermi level). Np/2 is the number of pairs of quasiparticles, taken in this work as the number of levels 
between the Fermi and the first level of the spectrum. 

For convenience, we have adopted the prescription of Ref. @, 0, which has been widely used for realistic potentials 
such as the Woods-Saxon potential (used here) or the folded- Yukawa potential Q . In this prescription, the force of 
the pairing is deduced from the empirical value of the gap A = 12/ VA and from N p (see text just above): 

i „^,„ (Jk. 



a *n^J (28) 

Here g(A) denotes the smoothed level density determined from the Strutinsky's procedure or by a semiclassical method 
as in the present work. The nonlinear system is solved by successive iterations until a given precision. At each iteration, 
we deduce the occupation probabilities from new couple A and A : 

' ' 1- ; -~ X | , fc = l, N p (29) 
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conversely, from the "new" occupations amplitudes (uk,Vk) we deduce the "new" gap: 



N P 

A = G^2u k vk avec u\ = 1 - v\ (30) 



and so on 

For one kind of particles, the pairing correction to the liquid drop model is defined as: 



pairing — corr 

were 



SPpairing-correct. (N Or Z, ft, 7) = P - P (31) 



P = Y? v ie k - - E 2efc ( 32 ) 
fe=i fe=i 

is the usual energie for a correlated system of fermions, and 

P = -i 9sc (A)A 2 (33) 

is its smooth part (i.e., without shell effects) assumed already contained in the liquid drop model 
Finally, with obvious notation, the potential energy of deformation can be summarized as follows: 

E def [N, Z- 0, 7) = E LD (N, Z; 0, 7) + SE SC {N; p, 7) + 5E SC (Z; /3, 7) 
+ 5P(N;P ll ) + 5P{Z ] p, 1 ) 

where the shell and pairing corrections are due to separates contributions of neutrons and protons. 
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VI. HANDLING AND NUMERICAL CHECKING OF THE ASSOCIATED FORTRAN CODE 



A. The non-deformed (spherical) case 

Two codes have been built for calculating the semiclassical energy. The first code is based on the general deformed 
case which consists of three fold integral (subroutine scdefor) and the second can only be used for the spherical shape 
with a one dimensional integral (subroutine sclspherl). Then, it is possible to make a cross checking in the spherical 
(non-deformed) case. To make further comparisons with other works, we have chosen the same examples as those of 
the Ref. [j| . The different contributions to the semiclassical energy Eq. (fl"0"l) are detailed in the following tables: 



TABLE I: Comparisons between our codes and other papers for the non-deformed case. The energies are given in 

MeV. 





N = 126; Vo = 44MeV; R v = Rso = 7.52/m; a v = a so = 0.67/ra; Hj = -0.7491; Z = arbitrary 


routine 


E sc 


A 


\N 


E% 


£°! 


E so 


K 


Ef° 


scdefor (present code) 
sclspherl (present code) 
Ref. [5] 


-2282.56 
-2282.61 
-2282.5 

TABLE 


-5.7476 
from scdefor 
-5.7323 

[I: Analog c; 


-724.21 
724.21 
-722.27 

ilculatioi 


1590.51 
1590.51 
1592.42 

is as in 


-50.90 
-50.87 
-50.87 

the tab 


24.08 
24.08 
23.92 

lefflfo 


from sclspherl 

-1.07 

-1.07 

r another exar 


from sclspherl 

-4.27 

-4.24 

uple 




N = 184; Vo = 43MeV; Rv = Rso = 8.48/m; a v = a so = 0.67/m; xj = -0.7321; Z = arbitrary 


routine 


Esc 


A 


XNo 


E% 


£°i 


E™ 


El 


Ef° 


scdefor (present code) 
sclspherl (present code) 
Ref. [5] 


-3228.76 
-3228.78 
-3230.0 


-4.9328 
from scdefor 
-4.9281 


-907.63 
-907.63 
-906.77 


2359.67 
2359.67 
2360.52 


-61.27 
-61.26 
-61.24 


30.12 
30.12 
30.13 


from sclspherl 

-2.62 

-1.47 


from sclspherl 

-4.76 

-4.97 



The numerical values of the parameters of the potential are displayed in the tables themselves. These calculations 
are performed for neutrons. The dependence on the proton number appears only through the parameters of the 
woods-saxon potential. Appart from numerical uncertainties due to different numerical approaches, the results are 
found very close. 

B. The deformed case 

To our knowledge, semiclassical calculations for the Hamiltonian such as the one considered in this paper do not 
exist in the literature. For this reason the only way to test the code in the deformed case is to compare the results 
with those of the Strutinsky type. However, it is well known that the latter method often gives results with some 
uncertainty. Consequently, as demonstrated in Ref. [2], in performing these tests, we must keep in mind that the 
Strutinsky calculations are only approximation of the semiclassical limit. In this respect, the smallness of the relative 
error gives a good idea on the quality of the results The essential point is to verify that the code runs properly. In 
fact the code has been checked extensively a longtime ago. As examples, we give two deformed cases in fig. (flj. The 
parameters are given in the readme4.pdf file. 

It is very clear that an approximative plateau exists in the region hu < 7 < l.bhuj represented by a circle. For 
the order p = we do not obtain any plateau. In the region of the plateau the relative error is less than lMeV per 
1300 - 1500MW in the both cases. 



VII. DATA, INPUT AND OUTPUT OF THE CODE 

This program has been designed on the Compac Visual Fortran version 6.6.0 (optimized settings). In 
fact, the structure of the code is somewhat complicated. So, it is no need to give too much de- 
tails. The essential point is to handle the basic input data and to be able to read the desired 
data from the output files. The fortran source code denoted by "enerdef.f can be downloaded from: 
http: / /made, voila.fr /index.php?m=c9ae77e8&a=7d397569&share=LNK80764b6393d92f388 
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FIG. 1: (Bottom) Semiclassical and Strutinsky calculations of the energy as function of the smearing parameter 7 
for different orders of the curvature correction. The semiclassical energy as well as the sum of single-particle energies 
are given by a straight line. The calculations are made for N = 54 (Bottom) and N=80 (Top). 
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A. Input data 

1. Parameters of the Woods-Saxon potential (file WSjparameters.dat) 

The microscopic model and the associated FORTRAN code is the same as the one of Ref.[| and 0. Therefore 
the parameters of the woods saxon potential are read from the file parameters.dat. renamed in the present work as 
ws_par ameters . dat . 

See pdf file readmel_woods saxon parameters 

2. Other input data (Beginning of the main program) 
See pdf file readme2_input data. 

They must be precised at the beginning of the main program: 
nmax=10 to 20 is linked to the size of the oscillator basis 
iuno=l (for single deformation) or (for lattice mesh points) 
if iuno=0 the three following data must be precised: 

betamax=0.0 up to about 1.0 is the maximal value of the parameter beta 
ibetapoints= is the number of points (minus one) in the beta direction 
igamapoints=is the number of points (minus one) in the gamma direction 

3. "Manual" input data (Keyboard) 

The kind of nucleons, the number of protons and the number of neutrons have to be entered manually on the 
keyboard. 

If iuno=l the deformation must also be precised in the terminal (do not forget that the deformation parameters 
are real quantities) 

4- Liquid drop data (Module liquid drop) 

The fissility parameter and other miscellaneous data for the liquid drop model are fixed in the subroutines eld, 
bbs,bbc in the module "liquid drop". 

5. Strutinsky calculations 

Additionally, this code is able to perform Strutinsky calculations. Two routines are devoted to this task. The first 
(Nstrutinsky) solves the Fermi level. The second (Strutinsky) calculates the smooth energy once the fermi level is 
known from the first routine. The essential points are the following for the rwo routines -(see readme3_strut.pdf file): 

ggam (input) = is the smearing parameter (ggam = 7 = Al.A~ 1 / 3 MeV) 

the numbers 0,8,16,18 (input, up to 18) = correspond to the curvature correction of the shell correction= does not 

exceed 18 (here four calculations are done). 

rnumbO, rnumb8, etc. ..(output for Nstrutinsky) = number of particle found after solving equation=checking 
hnew0,hnew8,.... (output for Nstrutinsky and input for Strutinsky)= Fermil level for different orders of the curv. 

correct. 

res0,res8,....= shell correction for different orders of the curv. correct. 

The code performs shell corrections in loop do for several values of ggam and four values of the order of the curvature 
correction. 



B. Data checkings 

In addition, the input and output data for the checkings are detailed in the readme4.pdf. file. 
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VIII. OUTPUT DATA 



The files eigenvalues and eigenvectors give the solutions of the Schrodinger equation. All results are given separately 
for neutrons and protons. Due to the coulomb interaction, the calculations in the proton case are significantly slower. 
However for a family of isotope the calculations for the protons must be taken only once. 

The files: 

del_n.dat, del_p.dat 
eldm_n.dat, eldm_p.da 
epot_n.dat, epot_p.dat 

give in the third column respectively the gap parameter, the energy of the liquide drop model and the deformation 
energy (all in MeV) for neutrons (_n) and protons (_p). The two first columns specify the deformation in the 
sextan, beta-gamma. Gamma is given in degrees. 

The files control results_n.dat and control results_p.dat give some details of the calculations. 

The File 2000n.dat (neutrons) or 2000p.dat (protons) gives the shell correction (columns 2 to 5). Each column 
corrresponds to a given order of the shell correction. Each row corresponds to a given value of the smearing parameter. 
The first column gives the smearing parameter (in hW units) and the last column gives the semiclassical value of the 
energy. 
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